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Abstract 

This paper presents a hybrid Godunov method for three-dimensional radiation hydro- 
dynamics. The multidimensional technique outlined in this paper is an extension of 
the one-dimensional method that was developed by Sekora & Stone 2009, 2010. The 
earlier one-dimensional technique was shown to preserve certain asymptotic limits and 
be uniformly well behaved from the photon free streaming (hyperbolic) limit through 
the weak equilibrium diffusion (parabolic) limit and to the strong equilibrium diffusion 
(hyperbolic) limit. This paper gives the algorithmic details for constructing a multidi- 
mensional method. A future paper will present numerical tests that demonstrate the 
robustness of the computational technique across a wide-range of parameter space. 



1 Radiation Hydrodynamics 

The purpose of this paper is to extend the ideas of Sekora & Stone 2010 to radiation hy- 
drodynamical problems in multiple dimensions. Therefore, this paper is a continuation of 
that earlier work, where the system of equations for radiation hydrodynamics were non- 
dimensionalized with respect to the material flow scale [TDl [XT] . This scaling gives two 
important parameters: C = c/aoo which measures relativistic effects and P = a r T^/ PooaL 
which measures how radiation affects material dynamics. Additionally, a r = 87r 5 /c 4 /15c 3 /i 3 = 
7.57 x 10~ 15 erg cm~ 3 K~ 4 is a radiation constant, is a reference material temperature in 
the absence of radiation, and is a reference material density in the absence of radiation. 
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The full system of equations for radiation hydrodynamics in the mixed frame that is correct 
to 0(1/C) is: 
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For the material quantities, p is density, m is momentum density, p = (7 — l)e is pressure, E 
is energy density, and T is temperature. For the radiation quantities, E r is energy density, 
F r is flux, P r is pressure, f is the variable tensor Eddington factor that is used to close 
the hierarchy of radiation transport moment equations, \ is a parameter, and n is a unit 
normal vector aligned with the radiative flux [6l[9]. In the free streaming limit (optically 
thin regime), \ — > 1 such that f — > n Cg) n. Yet, in the equilibrium diffusion limit (optically 
thick regime), x ~ * 1/3 such that f — > |/, where / is the identity matrix. In the above 
modified Mihalas-Klein source terms, a a is the absorption cross section, a s is the scattering 
cross section, and a t = cr a + a s is the total cross section [T014T2] . 



The numerical approach of the hybrid Godunov method rests on understanding the balance 
law form of the above system of equations, where: 
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The quasi-linear form of this system of hyperbolic balance laws is: 



4 



M D Sekora 



( 





1 




to 


to 














\ 






-(7- 


3)u 


-in - i)v 


— (7 — l)w 


(7-1) 





to 










-til) 


V 




u 


to 








to 










-TO 


w 




to 


u 








to 












h- in- 


-l)u 2 


— (7 — l)uv 


— (7 — 


■yu 





to 













to 




to 





to 





c 













to 




to 





to 


^■fxx 





















to 








&fyx 











V 










to 








Cf zx 








/ 



/ 










1 





















\ 




— uv 


V 




u 















to 










^V 2 - v 2 


-(7- 


l)u 


-(7- 


3)v 


"(7 


- l)w 


(7-1) 





to 










— vw 







w 






V 

















V 


(^V 2 H) 


-(7- 


\)uv 


h- in- 


-l)v 2 


-(7 


- l)vw 


7« 
















to 







to 


















c 







to 







to 

























to 







to 























V 


to 







to 












Cf zy 








/ 



/ 





to 





1 




to 














\ 




—uw 


w 





u 




to 





to 












—vw 


to 


w 


V 




to 





to 












Z^V 2 -w 2 


-(7- l)u 


-(7 - l)v 


-(7- 


3)w 


(7-1) 





to 












w (^v 2 -h) 


— (7 — l)uw 


— (7 — l)vw 




-l)w 2 


-fW 


















to 


to 


to 



















c 






to 


to 


to 










Cf xz 















to 


to 


to 










^■fyz 













V 


to 





to 










Cf zz 











/ 



Here, u = m x /p, v = m y /p, and w = m z /p are the velocities in the x, y, and z directions, 
respectively. V 2 = u 2 + v 2 + w 2 and H = ^ — ' 1 ^-V 2 is specific enthalpy. The Jacobians 
A x , A y , and A z have eigenvalues: X x = {0, u, u ± a, ±/i^ 2 C}, X y = {0, v , v ± a, i/^C}, 
and X z = {0,w,w ± a, ±/i/ 2 C}, respectively. However, one must account for how the stiff 
momentum and energy source terms affect the hyperbolic structure. 
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2 Overview of the Multidimensional Algorithm 

In radiation hydrodynamics, there are three important dynamical scales: the speed of sound 
(material flow), speed of light (radiation flow), and speed at which the source terms inter- 
act. Given such variation, one desires a numerical technique that treats the material flow 
explicitly, radiation flow implicitly, and source terms semi-implicitly. 



2.1 Effective CFL Condition 

Our primary interest when solving radiation hydrodynamical problems is higher-order reso- 
lution of material quantities while advancing the overall algorithm according to an effective 
CFL condition. This temporal advancement is associated with an effective CTU (corner 
transport upwind) scheme PQ: 

v 

At = p (It 

) \u(t,j,k)\+q cS (i,j,k) \v(%,j,k)\+q cfi (t,j,k) \w(i,j,k)\+q cf! (i,j,k) 
llicU H,i,fc I Ax ' Ay ' Az 

where At is the time step and v e [0, 1] is the CFL number. Ax = (x max — X nm>)/N^ ell , 
Ay = (y max - y m m)/N^ ell , and Az = (z max - z min )/N* eU are the spatial resolutions for 
a given number of computational grid cells in the x, y, and z directions, respectively. 
maxi tj}k {\u(i,j,k)\ + a eS (i,j,k), \v(i,j,k)\ +a e g(i,j,k), \w(i,j,k)\ + a cS (i,j,k)} is the maxi- 
mum material wave speed over all grid cells. Furthermore, a e s is an estimate of the radiation 
modified sound speed which is obtained by carrying out an effective eigen- analysis of the 
material Jacobian. This analysis is presented in a later section. One chooses a CTU-type 
time step over other alternatives (e.g., donor-cell time step) because of how one couples 
transport in the corners of the computational grid cells. For the duration of this paper, one 
assumes that the grid cells are cubic (Ax = Ay = Az). 

2.2 Algorithmic Steps 

After defining At, the algorithm loops over the following steps: 

1. Backward Euler Upwinding Scheme - implicitly advances the radiation quantities from 
time t n to time t n+ i\ 

(jpn pm pn pn \ , / Ptl+1 jpn+l pn+1 jpn+l\ 

\ r i r,xi r,yi r,zJ ' \ r ) r,x i r,y i r,z ) 



2. Modified Godunov Predictor Scheme - couples stiff source term effects to the hyperbolic 
structure of the balance laws for the material quantities and uses effective piecewise 
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linear extrapolation to spatially reconstruct material quantities at the left/right sides 
of cell interfaces in the x, y, and z directions {i ± 1/2, j ± 1/2, k ± 1/2}: 

i~i-m,n+l/2 fj-m,n+l/2 fjm,n+l/2 

U L/R,i+l/2,j,k> U L/R,i,j+l/2,k> U L/R,i,j,k+l/2 

3. Flux Function - evaluates the passage of material across cell interfaces using the 
left/right material states as well as an approximate Riemann solver to obtain: 

Fi+l/2,j,k = ^V^(U^i+i/2,j,k^R,i+l/2,j,k)) 
G^j+l/2,k = ^(^{^T,i,j+l/2,k^R,i,j+l/2,k)) 
-^ij',fc+l/2 = Hv^\U™ij,k+l/2i ^R,i,j,k+l/2)) 

4. CTU Correction - accounts for how the material quantities at the left/right sides of 
cell interfaces in one spatial direction are affected by the fluxes in the other two spatial 
directions: 

T~j-m,n+l/2 T7-m,n+l/2 
U L/R ~~ ^ U L/R 

5. Flux Function - evaluates the passage of material across cell interfaces using the cor- 
rected left/right material states and the algorithmic machinery of Step 3 above 

6. Modified Godunov Corrector Scheme - semi-implicitly advances the material quantities 
from time t n to time t n+1 : 

(p n ,m^m^m n z ,E n ) -> (p n+1 , m™ +1 , m™ +1 , m n z +1 , E n+l ) 

7. Apply boundary conditions 

8. Compute next time step At 

In the above expressions, U, U r , and U m represent all of the conserved quantities, radi- 
ation quantities (E r ,F rtX ,F rty ,F r}Z ), and material quantities (p,m x ,m y , m z ,E), respectively. 
Fi+i/2,j,ki Gij + \/2 t k, and Hij^+1/2 are fluxes directed across cell faces in the x, y, and z direc- 
tions, where 1Z represents the solution of a Riemann problem. The tilde that is above some of 
the left/right material states and fluxes designates quantities that have not yet been adjusted 
by the CTU correction. Cell centers are defined by three indices k), such that % ± 1/2, 
j ± 1/2, and k ± 1/2 represent the location of a cell interface to the right/left of i, j, and 
k, respectively, n designates the time discretization. Details about each step are explained 
in later sections. Lastly, the one-dimensional hybrid Godunov method of [18] was shown to 
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be consistent, stable, and accurate as well as coarsely gridded, well-balanced, and having 
some asymptotic preserving properties. For reasons cited in [18], these numerical properties 
should be able to be extended to the multidimensional algorithm. However, future tests will 
justify these claims. 

3 Backward Euler Upwinding Scheme 

This section presents the implicit scheme that advances the radiation quantities U r according 
to the material flow scale. The stability of explicit schemes (e.g., Godunov- type methods) is 
governed by the CFL condition which restricts the allowable time step according to the fastest 
characteristic speed. However, if a hyperbolic system consists of multiscale waves (e.g., 
radiation hydrodynamics where cja^ ~ 10 6 ), then explicit schemes can become inefficient. 
For these types of problems, implicit methods are useful. Following jU[7], one can construct 
implicit flux functions to approximate integrals at cell interfaces. In this context, a HLLE 
framework was implemented. 

3.1 HLLE Framework 

The HLLE scheme is based on estimating the minimum and maximum wave speeds (s m i n , s max ) 
These quantities are uniquely defined for each Riemann problem that is associated with each 
interface of a computational grid cell. The numerical flux in one direction is calculated using 
the following formula: 

F RLLE (n(U L , U R )) = i ((1 + C s ) (F{U L ) - s min U L ) + (1 - C s ) (F(U R ) - s max U R )) . (19) 

where C s = (s max + s min )/(s max — s min ). Defining the left/right states of the Riemann problem 
according to a first-order accurate (piecewise constant) reconstruction, forces the HLLE flux 
function to take the following form in each of the spatial directions: 

Fi+i/2,j,k(T^(U L ,i+i/2,j,k,U Rti+1 / 2 j,k)) ->■ Fi+yi,j,k('R'(Ui,j ) k, U i+1 j )k )), (20) 

Gij + l/2,k{7t{UL,i,j+l/2,k-,U R ^j + i/2,k)) — > G^ L +f/2,fc(^(^,j,fc> Uij + i : k)), (21) 

Hi,j,k+i/2(7t(UL,i,j,k+i/2,U Rt ij :k+ i/ 2 )) — »■ H it ^ k+ y 2 [JZ{JJi^ Uij } k+i))> (22) 

A first-order accurate, backward Euler-type algorithm was used because of total variation 
diminishing (TVD) conditions. These issues were explored by [HUE]. O ne makes the above 
explicit HLLE scheme implicit by defining the variables in the flux and source terms to be 



8 



M D Sekora 



at time t n+ \ such that the exact integral formulation of the conservative differencing is: 

U i,tk = U i,j,k ~ ^ ( F Sl/2, j,k(.^( U i J, L U i+l,j,k)) ~ F i^y%j,k( n { U i-l,j,h> U ij,k))) ( 23 ) 

^ ( rrHLLE (T?/TT n + 1 TT^ 1 \\ rrHLLE (<T?(TT n + 1 TT n + l \\\ 
~ \ H i,j,k+l/2\l<-\ U i,j,k > U i,j,k+l)) ~ H i,j,k-l/2\.'<-[ U i,j,k-V U i,j,k))) 
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3.2 Applying the Backward Euler HLLE Scheme 



If one considers only the radiation part of the equations for radiation hydrodynamics (Equa- 
tions H] and [5]) [HUE], then the variables, fluxes, and source terms are: 
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where the eigenvalues of the radiation subsystem in the free streaming limit (a a , cr t ~ 0{e)) 



are X x = {0,±/ix 2 C}, = {0,±/^ z C}, and X z = {0, ±fzl z C} for each of the spatial 
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directions. Given that the HLLE scheme uses minimum and maximum wave speeds to 
compute fluxes at cell interfaces, one defines the following equations: 



s x ■ 

"mm 



K,dhj, k) = -f„{ij, fc) 1/2 c, s * ax = a+h(< + i,j, k) = f„(i + fc) 1/2 c, 
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2 ' 



*Ln = KAh J,k) = -fyy(i,j,k)V 2 C, S max = A+ + l,k) = fyy(i,j + l,k)^C, 



/~is,y _ fyy{h 3 1) k) / fyy(i,j,k) I 
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r ^ /«(»,j,fc + l) 1/2 -/«(t,j,fc) 1/2 

y ' fc+1/2 /«(<,i,* + l) 1 / 2 + /„(i,i,*) 1 / 2 " 
Here, / w , and f zz arise from the closure relation P r = fE r and is either a user defined 
quantity or obtained by solving the radiation transport equation. If f varies spatially, then 
C s,x , C s,y , and C s ' z are non-zero. Defining or computing f(x, t) precedes the backward 
Euler update of U r . However, if f is assumed to be spatially and temporally constant, then 
C s,x ,C s ' y ,C S,Z = 0. Future work will update f(x, t) at each iteration by solving the radiation 
transport equation. 



3.3 Matrix Equation for the Radiation Components 

Inputting U r ' n+1 , F r (U r > n+1 ), G r {U r > n+1 ), H r (U r ' n+1 ) and S r (U m ' n ,U r > n+l ) into Equations 
12311261 gives the following four implicit difference equations: 
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+ d i ( x + C , S ,'/,fc+i/ 2 ) h^,3,k) - di (l - C 3 ;* fc _ 1/2 ) f yz (i,3,k) 
— - (m s + m x f X y(i, j, fe) + myfyy (i, j, fe) + m z f zy (i, j, fe)) + 



d 2 a a m i 



,3, fc) [l +di (l + C7f4f 1/2 ^ fc ) f xx (i,j, fc) 1/2 +di (l - C|l a j /2 ^_,.) f xx (i,j, fc) 1/2 

+ d! (l + C° : y +1/2 k ) f yy (l, j, fe)^ 2 +dl (l - C«;»_ 1/2ifc ) /„„(*, fe) 



j + 1 
3 + 1, 
1 (i, j, fe + 

1 (i,3,k + 



d i(i + c , s ;;, fc+ i/2) j. fc ) 1/2 + di (i - c , s ;;, fc -i/ 2 ) fc ) 1/2 + d 2- t ] 

fe) [dj (l -C'^^) fy X (i + l.j, fc)] 

fc) [-di(l-O i Ti/2.j,*)'~( < + 1 ^''*) 1/2 ] 
fc) [di (i - c";» +1/2ifc ) /»»(*, 3 + i, fc)] 

fc) [-dl (l - C e .'V +1/2 k ) f yy (i,3 + 1, fc) 1/2 ] 
1) [di {l-Cl : l k + 1/2 ) f yz (i,3,k + l)] 

1) [-di (l - C s ;J_ fc+1/2 ) /„(;, j, fc + l) 1 * 72 ] = F r ^(i, j, fc) + 



d 2 cr a m J/ T 4 



(31) 



pC 
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3,k- 
j, k - 
3 - 1 
j - 1 

J, fc) [ 



di 

di 

di 

di 

di 

di 

dl (l 

dl (l 

dl (l 
d 2 a"t 



i + c r;;, fc -i/ 2 )/"(».^ fc - i ) 1/2 ] 

1 + °iTl/!l,li) / *»f' , '- 1 ' fc )] 

1 + c u-vw) / »« (i ' j - 1 ' fc)1 1 

1 + ?-i/u») / '"( i - 1 ' 3 ''*) 1/J ] 
+ c ,Ti/2, 3 -, fe ) fc ) - d i (! - C r-l/2, 3 , fc ) /--(*. J- fe ) 

- C i*,'"+l/2,*) J. fe ) - dl (l - C i*,'"-l/2,*) f*V«,3,k) 

) fzz(i,3,k) - di (l - C 3 .'J k _ 1/2 ) f zz (i,j,k) 



»,j,fc+l/2 



- — - (m z + m x f X z(i, j, k) + m y f yz (i, j, fc) + m z f zz (i, j, fc)) H ' ° ' 1 

pC pC J 

j, fc) [l + di (l + C^j /2 /^(i, j, fc) 1/2 + di (l - C a .'* ll2 .^ f xx (i, j, fc) 1/2 

+ di (l + Of ■* +1/2 , )b ) /»»(i,J,fc) 1/2 +dl (l - C, S ;/_ 1/2 fc ) fyyii.j.k) 1 ' 2 

2 )f zz (i,j, fc) 1/2 + d 2 «r t ] 



E? +1 ( 
i+l 



l.J. 

1,3, 
J + 1, 
J + 1, 
J, fc + 



di(i + c s ;; ifc+1/2 )/„( l , J ,fc) 1 / 2 + di(i-c s 3ifc , , : 

di - - , , , , 



(i.j.fc + l) 



1 - C iTl/2,i,*)^( i + 1 '^ fc )] 



C 1 



t) S yv {i,3 + 1, *:) 



1/21 



d l ( 1 - C i, , /,* + l/2) / " (< ' : '' fe + 1) ] 



di (l 



i,j,fc + l/2 



)jWi,j,fe + l) 1/2 ] =F^(i,j,k). 



d 2 cr a m z T^ 
pC 



(32) 



where d\ = AtC/2Ax, d,2 = AtC All of the material quantities (p,m x ,m y ,m z , E) that 
appear in the above equations have the spatial index (i,j,k) and temporal index n. Ad- 
ditionally, the material temperature enters the above equations via the following relation: 

(E n ( JVf n \^ \ 

iff - \ > where M 2 = ml + ml + ml To under- 

stand how the multidimensional equations fit into the linear algebraic description Ax = b, 
it is important to understand that the algorithm cycles through the indices (i,j,k) in 
the following manner: f or(k=l ;N;k++){ f or (j=l ; N; j++) { f or (i=l ;N; i++) {...}}} . 
Furthermore, the algorithm cycles through the radiation quantities in the following order: 
(E r) F,^ X) F r , y) F r ^ z ). Assuming that N£ ell = N^ ell = N* ell , one expects the matrix A to have 
dimensions dim(A) = AN 3 x AN 3 . Moreover, the solution vector x will contain the following 
sequence of entries: 

x = [E? +1 (t,j, fc), F£\i,j, fc), F%\i,j, k), Fg\i,j, k), E^\ % + 1, j, k),..., 
E^\i + 2,j,k),..., Fg\NJ, fc), E? +1 (l,j + 1, fc), . . . , 
F^+ 1 (7V,7V,fc),^ +1 (l,l,fc + l),...] T . 
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If one assigns the character *asa placeholder for non-zero values, then one visualizes the 
structure of A which results from the above difference equations as a block banded matrix: 



0, *, 0, 0, *, <-> 0(iV 2 - N) <->,*, 0, *, 0, O 0(N - 
0, *, *, 0, 0, ++ 0(iV 2 - N) ++,*,*, 0, 0, ++ 0(iV - 
0, *, 0, *, 0, ++ 0(iV 2 - N) ++, *, 0, *, 0, ++ 0(iV - 
0, *, 0, 0, *, 0(iV 2 - N) ■(->, *, 0, 0, *, 0(iV - 



1) <->,*,*, 0, 0, *,*,*,*,*,*, 0, 0, ■<-> 0(iV - 
1) ++, *, *,0, 0, *, *, 0, 0, *,*, 0, 0, ++ 0(N - 
1) <+,*, 0, *, 0, *, 0, *, 0, *,0, *, 0, ++ 0(N - 
1) <->,*, 0, 0, *,*, 0, 0, *, *, 0, 0, *, ++ 0(JV - 



1) <->,*, 0, *, 0, O 0(N 2 - N) <->, *, 0, 0, *, 0, ... 
1) <->, *, *, 0, 0, ++ 0(N 2 - N) ++, *, *, 0, 0, 0, ... 
1) ++, *, 0, *, 0, ++ 0(N 2 - N) ++, *, 0, *, 0, 0, ... 
1) o, *, 0, 0, *, ++ 0(N 2 - N) «, *, 0, 0, *, 0, ... 



Since there are no nonlinearities in the radiation quantities for which root finding (e.g., New- 
ton's method) must be implemented, one casts these equations into a sparse matrix format 
that can be solved with iteration techniques from linear algebra such as the Jacobi, Gauss- 
Seidel, and multigrid methods as well as others because A exhibits diagonal dominance. 



4 Modified Godunov Predictor Scheme 

Given that the radiation quantities (E r ,F rjX ,F rj y,F rjZ ) are at time t n+ i, one computes the 
flux divergence ((V • F m ) n+1 / 2 ,(V • G m ) n+1 / 2 ,(V • H m ) n+1 ' 2 ) for the material quantities 
(p,m x ,m y ,m z , E) that are at time t n . Following the analysis of [T3 t[T7HT9] . one applies 
Duhamel's principle to the quasi-linear system of balance laws in Equations [bjlfTTl for only 
the material components. This technique defines the following system that locally includes 
in space and time the effects of the stiff source terms on the hyperbolic structure: 

DTI m f 8TI m 8TI m 8TI m \ 

"of = \~ A l^ ~ - + S m (U m >«, IT**)) , (33) 

where — ^- = ^ — h — h v ^ — h is the total derivative. L is a propagation 

operator that projects the dynamics of the stiff source terms onto the hyperbolic structure. 
A™ L = A™ — ul, A™ L = A™ — vl, and A™ L = A™ — wl are Jacobians associated with 
Lagrangian trajectories in the x, y, and z directions, respectively. Since the predictor scheme 
is a first-order accurate step in a second-order accurate predictor-corrector method, one 
chooses n = At/2 and the effective balance law becomes: 

dU m BU m BU m BU m 

~dT + A ^~dx~ + A ™^~dy~ + A ^~dz~ = J ( At /2)^ m (t/ m ' n , f/ r ' n+1 ), (34) 
where A™ eS = l(At/2)A™ L + ul, A™ eS = l(At/2)A™ L + vl, and A™ eS = l(At/2)A™ L + wl. 
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4.1 Applying the Modified Godunov Predictor Scheme 

If one considers only the material component in Equations [I]l3j then the variables, fluxes, 
and source terms are: 



U 



( 


p 


\ 




m x 






m y 






m z 




K 


E 


> 



H{U) 



( 



F{U) 



\ 



in 



m x 

-y+P 

m x m y 
P 

m x m z 



( 



G(U) 



m y 

m y m x 
„ P 

P V ~+P 

m y m z 



\ 



in; 



m z m x 
P 

'y 



m z m v 



-f+P 



S m (U) 



\ 





_ FS F X 

-FS p y 
-FS F * 
-FCS E 



\ 



In order to compute X, one first computes Vu m S m (U). Therefore: 



/ 

















\ 




-FS F * 














-FSp y 


-FS% X 


-FS% y 


-FSi\ 








-FS F * 


~^Sm x 






-P^E Z 




V 


-FCS F 


-FCS* x 




-PCS* 


-PCSf 


) 



V Um S m {U) 



where the partial derivatives are presented in Appendix 1. 



(35) 



4.1.1 Simplifying V Um S m (U) 

In its current form, Vu m S m (U) in Equation [351 leads to a propagation operator X that is 
difficult to work with algebraically. By inspection, the material momentum source terms are 
not the dominant factors defining the stiffness associated with the matter-radiation coupling, 
such that — PS F < 0(1) even in the strong equilibrium diffusion limit. Additionally, one 
finds that the derivative of the material momentum source term with respect to the variables 
U m has the following magnitude — PS^ m < 0(1). Therefore, — PS F can be included as a 
body force (e.g., gravity). 

It is the material energy source term — FCS E that defines the stiffness associated with the 
problem. By inspection of the contributing terms, — FCS E < 0(C) in the strong equilibrium 
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diffusion limit. Additionally, one finds that the derivative of the material energy source term 
with respect to the variables U m has the following magnitude — FCS§ m < 0(C 2 ). Therefore, 
one only needs to use — FCS E to define Vu™S m (U), such that: 



/ 



V Um S m {U) 










\ -FCS* 















-PCS£ 

















\ 



(36) 



-PCS* J 



Vu™S m (U) is further simplified by examining S^, S^ x , S^ y , S^ z , and S§ in the equilibrium 
diffusion limit and neglecting terms that have magnitudes of or less than (9(C). Therefore: 



'TO X 



Aa a T z ( 7 - 1) 
Aa a T 3 ( 7 - 1) 
4<7 a T 3 ( 7 -l) 
4<x a T 3 (7-1) 
4<7 a T 3 ( 7 -l) 



—E M< 



P z 

-m x 

P 2 
-my 

P 2 

-m z 



(37) 
(38) 
(39) 
(40) 
(41) 



It is important to note that these partial derivatives have the same stiff magnitude 4<r a T 3 (7 — 1). 
This insight simplifies algebraic manipulation. 



If V Um S m (U) is diagonalizable, then V Um S m {U) = RDR- 1 . Here, D = diag(0, 0, 0, 0, -PCS'!) 
and R is a matrix whose columns are the right eigenvectors. Below, one sees how the stiff 
magnitudes cancel out: 





-s§ 


qE 


_oE 
m y 


cE 










( 


qE 


qE 


qE 

m y 


qE 

rriz 






\ 







SZ 




SZ 



Sf 

1 






oE 
°E 




qE 

E 




qE 

°E 




qE 

E 
1 












1 




















1 










V 



1 


1 














1 


J 




V 




qE 

qE 
E 


1 

qE 
°m x 
qE 
E 




qE 

qE 
E 




qE 

qE 
E 




1 


/ 



(42) 
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4.1.2 Propagation Operator I 

Since one is considering a modified Godunov scheme with a predictor step of At/2 j!3j : 



At/2 
( 



At/2 



rV um S m (U) dr 



1 











o\ 





1 

















1 

















1 







(a-l)ff 

E 


(a-l)S 

E 


qE 

(«-!)% 





/ 1 







a) 



m2 

„2 




1 




- a) 



'1 





1 


- a) 







1 

7 p 



\ 




a y 



(43) 



(44) 



(45) 



where a = (l - exp(-PGSf At/2)) / (PCSfAt/2). Since Sf > across all relevant pa- 
rameter space, < a < 1. This property is important when considering stability and the 
subcharacteristic condition which is discussed later in this paper. 



4.2 Effective Material Jacobians 



Am Am Am 
^-x,effi ^Veff' ^eff 



The effects of the stiff source terms on the hyperbolic structure are accounted for by trans- 
forming to a moving-mesh (Lagrangian) frame (A™ L = A™ — ul, A 



y,L 



A™ - vl, A. 



z,L 



A m — wl), applying the propagation operator I to A™, and transforming back to an Eulerian 
frame, such that the effective material Jacobians (A™ eS = XA m L + ul, A™ eG = 1A™ L + vl, 
1-A™ L + wl) are given by [13J: 



Am 



-uv 

-UW 



-(7 - 3)t 



-(7 - l)u 2 + (^3_ + |1 



\ 

— (7 — l)v — (7 — (7 — 1) 

u 

u 

-(7 — l)uv —(7 — 7u j 



/ 1 

v u 

— (7 — l)u —(7 — 3)u — (7 — l)w (7 — 1) 

TJJ V 

I " [tV 2 - aH - (1 - a) (-^ + Iv 2 )] -(f-l)tM) -(7 " 1)« 2 + c*H + (1 - a) ( ^ + i^ 2 ) _( 7 _!)„«, 7 „ 



— liU 

izly 2 _ 
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— uw 

— vw 
L -V 2 -w 2 

-(1-")(^T + ^ 2 )] 



-(7 
-(7- 



- l)u 
l)uw 



~(7 
-(7- 



- l)v 

l)vw 



-(7 - l)w 2 + 



-(7 - 3)1/ 
aH + (1 - < 



0(^T + ^ 2 ) 







(7 - 1) 



These Jacobians have the following eigenvalues: A™ eff r + y = {u— a e s, u, u+a e s}, A™ cff r + , 
{v — a e ff, v,v + a c fi}, and A™ eff r + | = {w — a e ff, w,w + a e s}, respectively. Here, the effective 
sound speed a e g (i.e., the radiation modified sound speed) is: 

-.2 _ 1 ~ l-r/2 , „./ i\£V , n „.\ (m , J - ! TA 2 



< s = --^-V + a{ 7 - 1)H + (1 - a) \T + -L-^-V* J (46) 

= a — + (l-a)T (47) 
P 

= (a( 7 -l) + l)£, (48) 
P 

where T = p/p because of the choice of non-dimensionalization. Here, one notices that 
H, (T + (7 — l)V 2 /2) > across all relevant parameter space such that the effective sound 
speed a e ff admits the following limits: 

- FCS§ ^O^a^l a 2 cS -»■ -^—r-^V 2 + (7 - 1)# = — (adiabatic) (49) 

z p 

-FC^-oo=>„-0 «4^r = E (isothe^al). (50) 

When examining Equations H9] and [501 one sees that the subcharacteristic condition for 
material wave speeds is satisfied in each spatial direction, such that [13J: 



\ m 

A x,- 


= u - a ad < A™ effj _ 


= u 


— a c ff < A™ 


\ m 

— A x,eS,0 


— u < A™ cff + 


= u- 


- a cff < A™ eff!+ = u H 


- a a( j, 


(51) 


-V- 


- = u - a ad < A™ eff] _ 


= V 


- a ff < A™ 


\ m 

— A y,eH,0 




= v -+ 


- a eS < A™ cffi+ =dt 


Oad, 


(52) 




= tu - a a d < A™ eff _ 


= w 


— a c ff < A™ 


— A z,eff,0 


= ^ < A™ ff , + 


= w - 


f a cS < A™ cff + = w 


-1- a a d, 


(53) 



This condition is necessary for the stability of the system and guarantees that the numerical 
solution tends to the solution of the equilibrium equation as the relaxation time tends to 
zero. Additionally, the structure of the equations remains consistent with respect to clas- 
sical Godunov methods. Therefore, the CFL CTU time-stepping condition applies. Lastly, 
the right material eigenvectors BJP X y . eS (stored as columns) and left material eigenvectors 
^Tx y z\ eff (stored as rows) are given in Appendix 2. 
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4.3 Computing Left/Right States 

In the modified Godunov predictor scheme one uses effective piecewise linear extrapolation 
to spatially reconstruct material quantities at the left/right sides of cell interfaces. This 
technique is shown in the following relations for each spatial direction: 

^£v2,i,fc = u z% + (54) 



+ 





At 




Ax 






m 


i / 


At 



Vw&* = U ?^,k 2 V + ^ A *M U Z?,i,k)) PK^x^) (55) 

I f Qrn(Tjm,n Tpr,n+\ \ 

^£&. = u3i + \(i - ^yKA u ^)) p K AU S) (56) 

+ ^(f)^(^,^ 1 ), 

f T m,n+l/2 _ TT m,n ^ f j . A.t m (T jm,n \\ py ( \jjm,n \ /r 7 \ 

U R,i,j+l/2,k - U i,j+l,k 2 I ^A/.cfrl^U+l.fcJ J ^A^^iJ+l./J I ' J 



/-r,n+l \ 



+ f x (f) sm( ^'^ +1) ' 

Un£Sil/2 = U^ +1 - l(l + ^A^U^Sjn{^^.i) (59) 

i T ( 1 cm(Tjm,n jjr,n+l \ 

+ ~2 \~2~J l U i,j,fc+l' U i,j,k+V- 

where p^ ,v,z ^ is a slope limiting function used to eliminate spurious oscillations in the x, 
y, and z directions, rey, and?/, and?/, and z directions, rey, andy, and z directions, rey, 
and z directions, respectively. Slope limiting is performed for each of the material quan- 
tities. Although most slope limiters can be used, this algorithm employs the extremum- 
preserving [2l[T6] and traditional van Leer limiters (also referred to as the MUSCL limiter). 
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These techniques can be implemented either componentwise or across characteristic fields. 
After reconstructing the material quantities in each of the spatial directions, an approximate 
Riemann solver evaluates the passage of material into and out of each computational cell by 
using the material states that are to the left/right of the cell interfaces [T5]. It is important 
to emphasize that these flux functions do not directly account for the influence of radiation 
on the material quantities. Rather, the radiation effects are included via the source terms, 
propagation operator, and effective material Jacobian. 



5 Corner Transport Upwind Correction 

Before advancing the material quantities from time t n to time t n +i, one accounts for how 
the left/right states (U? /R>i+l/2>jtk , Ul/r^+w U™ /R .^ k+1/2 ) and thus the fluxes {F™ 1/w 
G™j+i/2k> Hi] k+1/2) are affected by transport in the other spatial directions. In particu- 
lar, one corrects for material propagating across the corners of a computational cell in the 
following manner pQ: 



TT m,n+l/2 _ jj-m,n+l/2 ^ L I ^ m _ ^ m \ /^ n \ 

U L,i+l/2,j,k - U L,i+l/2,j,k OA,, \ U i,j+y2,k W,i-l/2,fcJ l DU J 



At 

2A~y 
At 

2A~z 
At 



rrm rrm 

rL i,j,k+l/2 rL i,j,k-l/2 



jjm,n+l/2 _ jjm,n+l/2 / sim _ stm \ (n^ \ 

u R,i+i/2,j,k ~ u R,i+i/2,j,k 2Ay\ i+1 ' i+l/2 ' k ^t+ij-VWJ 1 



At 



2/± z \ 1A i+lj,k+l/2 J - t r+l,j,fc-l/2 



U L,i,j+l/2,k = U L,i,j+l/2,k 



U R,i,j+l/2,k = U R,i,j+l/2,k 



At 

2Ax 
At 



rpm rpm 
r i+l/2,j,k r i-l 



/2,j,k 



rrm TJ m 

2/\ z ^i,J,fc+l/2 n i,j,k-l/2 

( rpm _ rpm 

2~£^, \fi+l/2J+l,k r i-l/2,j+l,k 

At 



(62) 



(63) 



2Az 



rim _ rrm 

I1 i,j+l,k+l/2 n i,j+l,k-l/2 I 5 
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-1/2, k ; J 



TT m,n+l/2 _ Tf m,n+l/2 _ At / \ 

U L,i,i,fe+l/2 - U L,i,j,k+l/2 2Ax V i+1/2,j ' fc r i-l/2,i,fc y ^ 

2Ay V ^'.i-V 

T T m,n+l/2 _ f T m,n+l/2 _ At f ; _ v 

U R,i,j,k+l/2 ~ U R,i,j,k+l/2 2Ax V i + 1 / 2 .i>fe+l r i-l/2,i,fc+l J (uo ' 

At ( ^ m /im 

^i,j+l/2,fc+l ~ ^i,j-l/2,k+l 



With these corrected left/right material states, one again uses an approximate Riemann 
solver to evaluate the passage of material. 

6 Modified Godunov Corrector Scheme 

The time discretization for the source term is a single-step, second order accurate scheme 
based on [5j[T3l[Tl]. Given the material system of balance laws, one aims for a scheme that 
has an explicit approach for the flux divergence (V ■ F m + V ■ G m + V ■ H m ) and an implicit 
approach for the stiff source term S m (U). At each grid point, one solves the following 
collection of ordinary differential equations: 

^— = S m {U) - (V • F m ) n+1 ' 2 - (V ■ G m ) n+1/2 - (V • H m ) n+l ' 2 , (66) 

where the time-centered flux divergence terms are inputted from the predictor step and taken 
to be constant valued. Using Picard iteration and the method of deferred corrections, an 
initial guess for the solution to the collection of ordinary differential equations is: 

U = U m,n + At (/_ My Um S m {U)\ um , n)U r, n+ i)- 1 (67) 

( g m ^JJ rn ' n [7 r i n +l^ _ . prnyi+l/2 _ _ Qinyi+l/2 _ ^ _ jjm\n+l/2 j 



Here, Vu m S m (U) has the same functional form as that which was used to define the propa- 
gation operator X in a previous section. Therefore: 



(I - AtVu^S m (U)) 



\ AtFCS]? 




1 



AiPCS£ 






1 









1 

AiPCS^ 

Tfbz 










(68) 



1 + AtPCSf / 
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By inverting the above matrix, one finds: 



(/- AtVumSiU))- 1 



1 









1 









1 



-AtWCSB. 







1 

-AtPCS£_ 










(69) 



\ 1+AtPCSf 1+AtPCSf 1+AtPCSf 1+AtPCSf 1+AtPCSf / 

The error estimate e is the difference between the solution obtained from one iteration of the 
Picard technique where U is used as the starting value and the initial guess U : 

e(At) = U m ' n + ^(s m {U,U r ' n+1 ) + S m (U m ' n ,U r ' n+1 )} (70) 
- At ((V • F m ) n+1/2 + (V • G m ) n+1/2 + (V • ff m ) n+1/2 ) - U. 



Following Miniati & Colella 2007, the correction to the initial guess is given by: 

S(At) = (l- AtV um S m {U)\ (j ^ ur , n+i y 1 e(At). (71) 
Therefore, the material quantities at time t n +i are: 

U m ' n+1 = U + 6(At). (72) 



7 Conclusions and Future Work 

This paper presents the algorithmic details for constructing a hybrid Godunov method to 
solve three-dimensional radiation hydrodynamical problems. Careful consideration was taken 
when developing this technique such that one can compute numerical solutions for a host of 
physical phenomena (e.g., free streaming, weak equilibrium diffusion, and strong equilibrium 
diffusion limits). Additionally, the algorithmic ideas in this paper were cast in such a way 
so that they can be easily implemented in existing codes, particularly ones that carry out 
MHD calculations. Future papers will showcase (i) a series of multidimensional radiation 
hydrodynamical tests which demonstrates the robustness of the hybrid Godunov method 
and (ii) how to combine a numerical method for radiation hydrodynamics with a technique 
for updating the variable tensor Eddington factor f . 
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